\section{Normalization and Orthogonality}
Two modes, referenced by the indices $i=\{\ell_1,m_1,n_1\}$ and $j=\{\ell_2,m_2,n_2\}$ are defined to be \emph{orthonormal} if the following condition is true,
\begin{align}
\int_0^{z_0}\int_0^{2\pi}\vec{e}_{i}\times\vec{h}_{j}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}=\delta_{i,j}\label{eqn:defoforthognalitycyl}
\end{align}
where $\delta_{i,j}$ is the \emph{Kronecker delta} function and is defined as
\begin{align}
\delta_{i,j} = \left\{ \begin{array}{ll}
0 & \textrm{if $i\ne j$}\\
1 & \textrm{if $i=j$}
\end{array} \right.
\end{align}
Expanding the electric and magnetic basis vectors into $TE_z$ and $TM_z$ components yields the summation of four cases,
\begin{align}
\int_0^{z_0}\int_0^{2\pi}(\vec{e}_{Fi}+\vec{e}_{Ai})\times(\vec{h}_{Fj}+\vec{h}_{Aj})\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\label{eqn:defoforthognalitycylAF}
&=I_{FF}+I_{AA}+I_{FA}+I_{AF}
\end{align}
where
\begin{align}
I_{FF}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Fi}\times\vec{h}_{Fj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\label{eqn:IFFcyl}\\
I_{AA}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Ai}\times\vec{h}_{Aj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\label{eqn:IAAcyl}\\
I_{FA}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Ai}\times\vec{h}_{Fj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\label{eqn:IFAcyl}\\
I_{AF}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Fi}\times\vec{h}_{Aj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\label{eqn:IAFcyl}
\end{align}

Each of these four integrals will be analyzed separately.
\subsubsection{Case 1:  $TE_z\Leftrightarrow TE_z$}
Expanding the integrand and substituting (\ref{eqn:EFbasiscyl}) and (\ref{eqn:HFbasiscyl}) into (\ref{eqn:IFFcyl}) gives
\begin{multline}\label{eqn:IFFcyl2}
I_{FF}=\int_0^{z_0}\int_0^{2\pi}e_{F\phi i}h_{Fzj}\;({\rho}{\;d\phi}){\;dz}=\\ C_{f\ell_1,m_1,n_1}C_{f\ell_2,m_2,n_2}k_{f\rho{n_1}}k_{f\rho{n_2}}\sqrt{\frac{Z_{f{n_1}}}{Z_{f{n_2}}}}{\rho}I_{\Phi FF}I_{ZFF}
\end{multline}
where,
\begin{align}
I_{\Phi FF} &= \int_0^{2\pi}\Phi_F(\phi)_{\ell_1,m_1}\Phi_F(\phi)_{\ell_2,m_2}{\;d\phi}\label{eqn:IPhicyl}\\
&=\int_0^{2\pi}\cos\left(m_1\phi-\ell_1\frac{\pi}{2}\right)\cos\left(m_2\phi-\ell_2\frac{\pi}{2}\right){\;d\phi}
\end{align}
\begin{align}
I_{ZFF} &= \int_0^{z_0}Z_F(z)_{n_1}Z_F(z)_{n_2}{\;dz}\label{eqn:IZcylTE}\\
&=\int_0^{z_0}\sin\left(\frac{n_1\pi}{z_0}z\right)\sin\left(\frac{n_2\pi}{z_0}z\right){\;dz}
\end{align}
When $m_1\ne m_2$ and for any combination of $\ell_1$ and $\ell_2$,
\begin{align}
I_{\Phi FF}&=0 \label{eqn:IPHIcylequal0}
\end{align}
When $n_1\ne n_2$ then,
\begin{align}
I_{ZFF}&=0
\end{align}
which proves that the modes are orthogonal. In order to make the modes orthonormal, let $i=j$ so that (\ref{eqn:IFFcyl2}) reduces to,
\begin{align}
C_{f\ell,m,n}^2k_{f\rho{n}}^2{\rho}I_{\Phi FF}I_{ZFF}=1\label{eqn:exhcyltemp}
\end{align}
where
\begin{align}
I_{\Phi FF}&=\int_0^{2\pi}\Phi_F^2(\phi)_{\ell,m}{\;d\phi}\\
&=\int_0^{2\pi}\cos^2\left(m\phi-\ell\frac{\pi}{2}\right){\;d\phi}\\
&=\pi\left[1+(-1)^\ell\delta_{m,0}\right]\label{eqn:intPhisqr}
\end{align}
and where
\begin{align}
I_{ZFF}&=\int_0^{z_0}Z_F^2(z)_{n}{\;dz}\\
&=\int_0^{z_0}\sin^2\left(\frac{n\pi}{z_0}z\right){\;dz}\\
&=\left(1-\delta_{n,0}\right)\frac{z_0}{2}\label{eqn:intZsqr}
\end{align}
Finally solving for $C_f$ in (\ref{eqn:exhcyltemp}) and substituting (\ref{eqn:intPhisqr}) and (\ref{eqn:intZsqr}) yields,
\begin{align}
C_{f\ell,m,n} = \left[{\rho}{z_0}k_{f\rho{n}}^2\frac{\pi}{2}\left[1+(-1)^\ell\delta_{m,0}\right]\left[1-\delta_{n,0}\right]\right]^{-\frac{1}{2}}
\end{align}

\subsubsection{Case 2:  $TM_z\Leftrightarrow TM_z$}
\begin{align}
I_{AA}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Ai}\times\vec{h}_{Aj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\\&=
\int_0^{z_0}\int_0^{2\pi}-e_{Azi}h_{A\phi j}\;({\rho}{\;d\phi}){\;dz}\label{eqn:IAA}\\
&=C_{a\ell_1,m_1,n_1}C_{a\ell_2,m_2,n_2}k_{a\rho{n_1}}k_{a\rho{n_2}}\sqrt{\frac{Z_{a{n_1}}}{Z_{a{n_2}}}}{\rho}I_{\Phi AA}I_{ZAA}
\end{align}
Expanding the integrand and substituting (\ref{eqn:EAbasiscyl}) and (\ref{eqn:HAbasiscyl}) into (\ref{eqn:IAAcyl}) gives
\begin{multline}
\int_0^{z_0}\int_0^{2\pi}-e_{Azi}h_{A\phi j}\;({\rho}{\;d\phi}){\;dz}\\=  C_{a\ell_1,m_1,n_1}C_{a\ell_2,m_2,n_2}k_{a\rho{n_1}}k_{a\rho{n_2}}\sqrt{\frac{Z_{a{n_1}}}{Z_{a{n_2}}}}{\rho}I_{\Phi AA}I_{ZAA}
\end{multline}
where $I_{\Phi AA}$ is given in (\ref{eqn:IPhicyl})
\begin{align}
I_{\Phi AA} &= \int_0^{2\pi}\Phi_A(\phi)_{\ell_1,m_1}\Phi_A(\phi)_{\ell_2,m_2}{\;d\phi}\label{eqn:IPhicylA}\\
&=\int_0^{2\pi}\cos\left(m_1\phi-\ell_1\frac{\pi}{2}\right)\cos\left(m_2\phi-\ell_2\frac{\pi}{2}\right){\;d\phi}
\end{align}
and,
\begin{align}
I_{ZAA} &= \int_0^{z_0}Z_A(z)_{n_1}Z_A(z)_{n_2}{\;dz}\label{eqn:IZcylTM}\\
&=\int_0^{z_0}\cos\left(\frac{n_1\pi}{z_0}z\right)\cos\left(\frac{n_2\pi}{z_0}z\right){\;dz}
\end{align}
When $m_1\ne m_2$ and for any combination of $\ell_1$ and $\ell_2$ then $I_{\Phi AA}=0$. When $n_1\ne n_2$ then $I_{ZAA}=0$
which proves that the modes are orthogonal. In order to make the modes orthonormal, let $i=j$ so that (\ref{eqn:IAAcyl}) reduces to,
\begin{align}
C_{a\ell,m,n}^2k_{a\rho{n}}^2{\rho}I_{\Phi AA}I_{ZAA}
=1\label{eqn:exhcyltemp2}
\end{align}
where $I_{\Phi AA}$ reduces to,
\begin{align}
I_{\Phi AA} &= \int_0^{2\pi}\Phi^2_A(\phi)_{\ell,m}{\;d\phi}\\
&=\int_0^{2\pi}\cos^2\left(m\phi-\ell\frac{\pi}{2}\right){\;d\phi}\\
&=\pi\left[1+(-1)^\ell\delta_{m,0}\right]\label{eqn:intPhisqrA}
\end{align}
and 
\begin{align}
I_{ZAA}&=\int_0^{z_0}Z^2_A(z)_{n}{\;dz}\\
&=\int_0^{z_0}\cos^2\left(\frac{n\pi}{z_0}z\right){\;dz}\\
&=\left(1+\delta_{n,0}\right)\frac{z_0}{2}\label{eqn:intZsqrTM}
\end{align}
Finally solving for $C_a$ in (\ref{eqn:exhcyltemp2}) and substituting (\ref{eqn:intPhisqrA}) and (\ref{eqn:intZsqrTM}) yields,
\begin{align}
C_{a\ell,m,n} = \left[{\rho}{z_0}k_{a\rho{n}}^2\frac{\pi}{2}\left[1+(-1)^\ell\delta_{m,0}\right]\left[1+\delta_{n,0}\right]\right]^{-\frac{1}{2}}
\end{align}

\subsubsection{Case 3:  $TE_z\Leftrightarrow TM_z$}
Expanding the cross product of the integrand in (\ref{eqn:IFAcyl}) yields,
\begin{align}
I_{FA}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Ai}\times\vec{h}_{Fj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}=0\label{eqn:IAF}
\end{align}

\subsubsection{Case 4:  $TM_z\Leftrightarrow TE_z$}
\begin{align}
I_{AF}&=\int_0^{z_0}\int_0^{2\pi}\vec{e}_{Fi}\times\vec{h}_{Aj}\cdot\hat{\rho}\;({\rho}{\;d\phi}){\;dz}\\
&=\int_0^{z_0}\int_0^{2\pi}(e_{A\phi i}h_{Fzj}-e_{Azi}h_{F\phi j})\;({\rho}{\;d\phi}){\;dz}\label{eqn:IFA}
\end{align}
Expanding the integrand and substituting (\ref{eqn:EFbasiscyl}), (\ref{eqn:HFbasiscyl}), (\ref{eqn:EAbasiscyl}) and (\ref{eqn:HAbasiscyl}) yields,
\begin{align}
I_{AF}&=C_{a}C_{f}\frac{\sqrt{Z_a}}{\sqrt{Z_f}}\left[\frac{k_{f\rho}}{k_{a\rho}}I_{\Phi FA}I_{ZFA}-\frac{k_{a\rho}}{k_{f\rho}}I_{\Phi AF}I_{ZAF}\right]
\end{align}
where
\begin{align}
I_{\Phi FA}&=\int_0^{2\pi}\Phi_F(\phi)\Phi'_A(\phi)\;d\phi\\
&=-m_2\int_0^{2\pi}\cos\left(m_1\phi-\ell_1\frac{\pi}{2}\right)\sin\left(m_2\phi-\ell_2\frac{\pi}{2}\right)\;d\phi\\
I_{ZFA}&=\int_0^{z_0}Z_F(z)Z'_A(z)\;dz\\
&=-\frac{n_2\pi}{z_0}\int_0^{z_0}\sin\left(\frac{n_1\pi}{z_0}z\right)\sin\left(\frac{n_2\pi}{z_0}z\right)\;dz\\
I_{\Phi AF}&=\int_0^{2\pi}\Phi_A(\phi)\Phi'_F(\phi)\;d\phi\\
&=-m_1\int_0^{2\pi}\cos\left(m_2\phi-\ell_2\frac{\pi}{2}\right)\sin\left(m_1\phi-\ell_1\frac{\pi}{2}\right)\;d\phi\\
I_{ZAF}&=\int_0^{z_0}Z_A(z)Z'_F(z)\;dz\\
&=\frac{n_1\pi}{z_0}\int_0^{z_0}\cos\left(\frac{n_2\pi}{z_0}z\right)\cos\left(\frac{n_1\pi}{z_0}z\right)\;dz
\end{align}

\subsection{Power $TE_z$}
The total power is given by the \emph{Poynting theorem} as,
\begin{align}
P &= \frac{1}{2}\Re\left\{\iint_S\vec{E}\times\vec{H}^*\cdot d\vec{s}\right\}\\
&= \frac{1}{2}\Re\left\{\int_0^{2\pi}\int_0^{z_0}\vec{E}_\bot\times\vec{H}_\bot^*\; (\rho d\phi)dz\right\}
\end{align}
where $^*$ is the complex conjugate operator. 
Substituting (\ref{eqn:EFperpcyl}) and (\ref{eqn:EFperpcyl}) yields,
\begin{align}
P &= \frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R'_{\ell,m,n}(\rho)R_{\ell,m,n}^*(\rho)}{k_{\rho n}}\int_0^{2\pi}\int_0^{z_0}\vec{e}_{\ell,m,n}(\phi,z)\times\vec{h}^*_{\ell,m,n}(\phi,z)\; (\rho d\phi)dz\right\}\\
&=\frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R'_{\ell,m,n}(\rho)R_{\ell,m,n}^*(\rho)}{k_{\rho n}}\int_0^{2\pi}\int_0^{z_0}{e}_{\phi\ell,m,n}(\phi,z){h}^*_{z\ell,m,n}(\phi,z)\;(\rho d\phi)dz\right\}
\end{align}
Substituting the basis vectors (\ref{eqn:EFbasiscyl}) and (\ref{eqn:HFbasiscyl}) yields,
\begin{align}
P &= \frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R'_{\ell,m,n}(\rho)R_{\ell,m,n}^*(\rho)}{k_{\rho n}}C_{f\ell,m,n}C^*_{f\ell,m,n}k_{\rho{n}}k^*_{\rho{n}}\frac{\sqrt{Z_{fn}}}{\sqrt{Z^*_{fn}}}\rho I_{\phi\phi}I_{zz}\right\}\label{eqn:Pwrcyl}
\end{align}
where $I_{\phi\phi}$ and $I_{zz}$ are defined in (\ref{eqn:intPhisqr}) and (\ref{eqn:intZsqr}). From (\ref{eqn:exhcyltemp}) it is found that,
\begin{align}
I_{\phi\phi}I_{zz}=\frac{1}{C_{f\ell,m,n}^2k_{\rho{n}}^2{\rho}}
\end{align}
and then substituting,
\begin{align}
P=\frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R'_{\ell,m,n}(\rho)R_{\ell,m,n}^*(\rho)}{k_{\rho{n}}}\frac{C^*_{f\ell,m,n}}{C_{f\ell,m,n}}\frac{k^*_{\rho{n}}}{k_{\rho{n}}}\frac{\sqrt{Z_{fn}}}{\sqrt{Z^*_{fn}}}\right\}\label{eqn:Pwrcyl2}
\end{align}
Since,
\begin{align}
\frac{C^*_{f\ell,m,n}}{C_{f\ell,m,n}}\frac{k^*_{\rho{n}}}{k_{\rho{n}}}=1
\end{align}
then the total power power is given by
\begin{align}
P&=\frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R'_{\ell,m,n}(\rho)R_{\ell,m,n}^*(\rho)}{k_{\rho{n}}}e^{\angle{Z_{fn}}}\right\}\label{eqn:Pwrcyl3}\\
\end{align}
where,
\begin{align}
e^{\angle{Z_{fn}}}=\frac{\sqrt{Z_{fn}}}{\sqrt{Z^*_{fn}}}
\end{align}
It can be seen that the the total power is a summation of the power of each individual mode.

\subsection{Power $TM_z$}
%The total power is given by the \emph{Poynting theorem} in (\ref{eqn:PoyntingThereom}) and (\ref{eqn:PoyntingVector}).
Substituting (\ref{eqn:EAperpcyl}) and (\ref{eqn:EAperpcyl}) yields,
\begin{align}
\vec{S}(\phi,z) &= \frac{R_{\ell,m,n}(\rho)R'{_{\ell,m,n}^*}(\rho)}{k^*_{\rho n}}\vec{s}(\phi,z)\label{eqn:PoyntingVectorTM2}
\end{align}
where,
\begin{align}
\vec{s}(\phi,z) = \frac{1}{2}\vec{e}_{\ell,m,n}(\phi,z)\times\vec{h}^*_{\ell,m,n}(\phi,z)=\frac{1}{2}{e}_{\phi\ell,m,n}(\phi,z){h}^*_{z\ell,m,n}(\phi,z)\label{eqn:PoyntingVectorBasisTM}
\end{align}
Substituting the basis vectors (\ref{eqn:EAbasiscyl}) and (\ref{eqn:HAbasiscyl}) yields,
\begin{align}
P &= \frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R_{\ell,m,n}(\rho)R'{_{\ell,m,n}^*(\rho)}}{k^*_{\rho n}}C_{a\ell,m,n}C^*_{a\ell,m,n}k_{\rho{n}}k^*_{\rho{n}}\frac{\sqrt{Z_{an}}}{\sqrt{Z^*_{an}}}\rho I_{\phi\phi}I_{zz}\right\}\label{eqn:PwrcylTM}
\end{align}
where $I_{\phi\phi}$ and $I_{zz}$ are defined in (\ref{eqn:intPhisqr}) and (\ref{eqn:intZsqrTM}). From (\ref{eqn:exhcyltemp2}) it is found that,
\begin{align}
I_{\phi\phi}I_{zz}=\frac{1}{C_{a\ell,m,n}^2k_{\rho{n}}^2{\rho}}
\end{align}
and then substituting,
\begin{align}
P=\frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R_{\ell,m,n}(\rho)R'{_{\ell,m,n}^*(\rho)}}{k^*_{\rho{n}}}\frac{C^*_{a\ell,m,n}}{C_{a\ell,m,n}}\frac{k^*_{\rho{n}}}{k_{\rho{n}}}\frac{\sqrt{Z_{an}}}{\sqrt{Z^*_{an}}}\right\}\label{eqn:Pwrcyl2TM}
\end{align}
Since,
\begin{align}
\frac{C^*_{a\ell,m,n}}{C_{a\ell,m,n}}\frac{k^*_{\rho{n}}}{k_{\rho{n}}}=1
\end{align}
then the total power power is given by
\begin{align}
P&=\frac{1}{2}\Re\left\{\sum_{\ell,m,n}\frac{R_{\ell,m,n}(\rho)R'{_{\ell,m,n}^*(\rho)}}{k^*_{\rho{n}}}e^{\angle{Z_{an}}}\right\}\label{eqn:Pwrcyl3TM}\\
\end{align}
where,
\begin{align}
e^{\angle{Z_{an}}}=\frac{\sqrt{Z_{an}}}{\sqrt{Z^*_{an}}}
\end{align}
It can be seen that the the total power is a summation of the power of each individual mode.